Model for joint multispecies distributions

Using real data. For more infromation on how these data were obtain go to: http://panthera:8888/notebooks/external_plugins/paper3code/%5BD.E%5DCell%20Extractor.ipynb

Change to the working directory

Load libraries

The following libraries are needed for fitting the model

These libraries are used for visualising the posterior sample

Data exploration

The following line gives the structure of the dataframe.

Plotting the spatial distribution of the data set (data)

Here we show the spatial observations of Disocactus (id = 1) on the level = 1

Display a section (100:100) of the adjacency matrix

Fitting the model

Define the formulas for the ecological process ($P$) and the sampling effort ($S$). In this case the ecological process would be defined as:

$$ logit(P_y) = Elevation + Precipitation + MeanTemperature + G$$$$ logit(S) = Ecoregions + Typology + G$$

Note:

Elevation, Precipitation and Mean Temperature will be standardized (i.e. $\frac{X}{\sigma²} - \mu_X$ X $\in$ {Elevation, Precipitation, Mean Temperature})

wwf_mhtnam is the name of the column corresponding to the Terrestrial Ecoregions, tipolgia14 is the name of the column corresponding to the settlement typology.

Sample the posterior distribution generated by the model

We will sample 1000 iterations using four chains (no thinning) using Hamiltonian Monte Carlo.

Analysis of results

Calculate the summary statistics of the fixed effects (parameters beta), the mixed component between the ecological process $P_i$ and the sampling effort $S$, the spatial autocorrelation (alpha_car) and the overall variance (tau).

Data visualisation

Uncomment this to save the STAN posterior sample object

Saving the model and traces

import pickle _file = "/outputs/presence_only_models/multispecies_model.pkl" with open(_file,"wb") as f: %time pickle.dump([multispecies_model,fit3],f,protocol=-1)

Continue sampling after burnin

Test to unpickle an object

print(fit_continue.stansummary(pars='beta')) print(fit_continue.stansummary(pars='alpha')) print(fit_continue.stansummary(pars='alpha_car')) print(fit_continue.stansummary(pars='tau'))